%% Sup figure: Turning rate ramp
load('./SimData/Parameters.mat')
load('./SimData/MatchedValues.mat')
load('./SimData/TurningRateRamp.mat')
temp=fields(Parameters);
for i = 1:length(temp)
    eval([temp{i} '= Parameters.' temp{i} ';'])
end
row = find(V_r(:,1) == Vr_match(5));

cmap = cool(256);
figure; hold on; box on;
colormap cool

log_turning = log10(TurningRates);
for ii = 1:length(TurningRates)
    
    temp = log_turning(ii);temp = (temp + abs(min(log_turning)))./(max(log_turning)-min(log_turning));
    
    plot(bin_centers, flip(TurningRateRamp(ii,:)), 'color', cmap(floor(255*temp)+1,:))
    
    if ii == length(TurningRates)
        V = (1+bin_centers*(V_r(row,ii)-100)/100);
        k = 1./trapz(bin_centers, 1./V);
        rho_s = k./V;
        plot(bin_centers, flip(rho_s), '--k', 'linew', 2)
        writematrix(flip(rho_s),'ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', ['C' num2str(ii+3)])
        writematrix("0_theoretical",'ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', ['A' num2str(ii+3)])
    end
    writematrix(flip(TurningRateRamp(ii,:)),'ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', ['C' num2str(ii+2)])
end
cbar = colorbar;
caxis([min(log_turning) max(log_turning)])
cbar.Ticks = [-3 -2 -1];
xticks([0 0.25 0.5 0.75 1])
xlabel('position')
ylabel('Density (PDF)')
title(cbar, 'log10 Turning rate (1/\tau)')
set(gca,'linewidth', 1.5)
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [6 4]);

writematrix('position (um)','ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', 'B1')
writematrix(bin_centers*1000,'ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', 'C1')
writematrix('turning rate (s^-1)','ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', 'A2')
writematrix(1./list_Tau_visc','ExtDataFigure3.xlsx', 'Sheet', 'Ext3', 'Range', 'A3')
print('./TurningRateRamp/phasematrixverticallinecut', '-dpdf', '-r1200', '-bestfit')

%% FUNCTION

function p = FPdistribution(theta, t, Dr, tau)
Nt = length(t);
Ntheta = length(theta);
dt = abs(t(2)-t(1));
dtheta = abs(theta(2)-theta(1));
p0 = 1/2/pi*ones(size(theta));

p = zeros(Nt,Ntheta);

for ii = 1:Nt
    
    if ii == 1
        p(ii,:) = p0;
    else
        p(ii,2:end-1) =  p(ii-1,2:end-1) ...
            + dt*(   Dr * (p(ii-1,3:end) - 2*p(ii-1,2:end-1) + p(ii-1,1:end-2))/dtheta^2 ...
            + (1/tau) * (p(ii-1,3:end).*sin(theta(3:end)) - p(ii-1,1:end-2).*sin(theta(1:end-2)))/2/dtheta);
        
        p(ii,1) =  p(ii-1,1) ...
            + dt*(   Dr * (p(ii-1,2) - 2*p(ii-1,1) + p(ii-1,end))/dtheta^2 ...
            + (1/tau) * (p(ii-1,2).*sin(theta(2)) - p(ii-1,end).*sin(theta(end)))/2/dtheta);
        p(ii,end) =  p(ii-1,end) ...
            + dt*(   Dr * (p(ii-1,1) - 2*p(ii-1,end) + p(ii-1,end-1))/dtheta^2 ...
            + (1/tau) * (p(ii-1,1).*sin(theta(1)) - p(ii-1,end-1).*sin(theta(end-1)))/2/dtheta);
    end
    
    
end
end
%% FUNCTION
function [param,stat]=sigm_fit(x,y,fixed_params,initial_params,plot_flag)
% Optimization of parameters of the sigmoid function
%
% Syntax:
%       [param]=sigm_fit(x,y)
%
%       that is the same that
%       [param]=sigm_fit(x,y,[],[],[])     % no fixed_params, automatic initial_params
%
%       [param]=sigm_fit(x,y,fixed_params)        % automatic initial_params
%       [param]=sigm_fit(x,y,[],initial_params)   % use it when the estimation is poor
%       [param]=sigm_fit(x,y,fixed_params,initial_params,plot_flag)
%
% param = [min, max, x50, slope]
%
% if fixed_params=[NaN, NaN , NaN , NaN]        % or fixed_params=[]
% optimization of "min", "max", "x50" and "slope" (default)
%
% if fixed_params=[0, 1 , NaN , NaN]
% optimization of x50 and slope of a sigmoid of ranging from 0 to 1
%
%
% Additional information in the second output, STAT
% [param,stat]=sigm_fit(x,y,fixed_params,initial_params,plot_flag)
%
%
% Example:
% %% generate data vectors (x and y)
% fsigm = @(param,xval) param(1)+(param(2)-param(1))./(1+10.^((param(3)-xval)*param(4)))
% param=[0 1 5 1];  % "min", "max", "x50", "slope"
% x=0:0.1:10;
% y=fsigm(param,x) + 0.1*randn(size(x));
%
% %% standard parameter estimation
% [estimated_params]=sigm_fit(x,y)
%
% %% parameter estimation with forced 0.5 fixed min
% [estimated_params]=sigm_fit(x,y,[0.5 NaN NaN NaN])
%
% %% parameter estimation without plotting
% [estimated_params]=sigm_fit(x,y,[],[],0)
%
%
% Doubts, bugs: rpavao@gmail.com
% Downloaded from http://www.mathworks.com/matlabcentral/fileexchange/42641-sigmoid-logistic-curve-fit

% warning off

x=x(:);
y=y(:);

if nargin<=1 %fail
    fprintf('');
    help sigm_fit
    return
end

automatic_initial_params=[quantile(y,0.05) quantile(y,0.95) NaN 1];
if sum(y==quantile(y,0.5))==0
    temp=x(y==quantile(y(2:end),0.5));
else
    temp=x(y==quantile(y,0.5));
end
automatic_initial_params(3)=temp(1);

if nargin==2 %simplest valid input
    fixed_params=NaN(1,4);
    initial_params=automatic_initial_params;
    plot_flag=1;
end
if nargin==3
    initial_params=automatic_initial_params;
    plot_flag=1;
end
if nargin==4
    plot_flag=1;
end

if exist('fixed_params','var')
    if isempty(fixed_params)
        fixed_params=NaN(1,4);
    end
end
if exist('initial_params','var')
    if isempty(initial_params)
        initial_params=automatic_initial_params;
    end
end
if exist('plot_flag','var')
    if isempty(plot_flag)
        plot_flag=1;
    end
end

%p(1)=min; p(2)=max-min; p(3)=x50; p(4)=slope como em Y=Bottom + (Top-Bottom)/(1+10^((LogEC50-X)*HillSlope))
%f = @(p,x) p(1) + (p(2)-p(1)) ./ (1 + 10.^((p(3)-x)*p(4)));

f_str='f = @(param,xval)';
free_param_count=0;
bool_vec=NaN(1,4);
for i=1:4;
    if isnan(fixed_params(i))
        free_param_count=free_param_count+1;
        f_str=[f_str ' param(' num2str(free_param_count) ')'];
        bool_vec(i)=1;
    else
        f_str=[f_str ' ' num2str(fixed_params(i))];
        bool_vec(i)=0;
    end
    if i==1; f_str=[f_str ' + (']; end
    if i==2;
        if isnan(fixed_params(1))
            f_str=[f_str '-param(1) )./ (   1 + 10.^( ('];
        else
            f_str=[f_str '-' num2str(fixed_params(1)) ')./ (1 + 10.^(('];
        end
    end
    if i==3; f_str=[f_str ' - xval ) *']; end
    if i==4; f_str=[f_str ' )   );']; end
end

eval(f_str)

[BETA,RESID,J,COVB,MSE] = nlinfit(x,y,f,initial_params(bool_vec==1));
stat.param=BETA';

% confidence interval of the parameters
stat.paramCI = nlparci(BETA,RESID,'Jacobian',J);

% confidence interval of the estimation
[stat.ypred,delta] = nlpredci(f,x,BETA,RESID,'Covar',COVB);
stat.ypredlowerCI = stat.ypred - delta;
stat.ypredupperCI = stat.ypred + delta;

% plot(x,y,'ko') % observed data
% hold on
% plot(x,ypred,'k','LineWidth',2)
% plot(x,[lower,upper],'r--','LineWidth',1.5)

free_param_count=0;
for i=1:4;
    if isnan(fixed_params(i))
        free_param_count=free_param_count+1;
        param(i)=BETA(free_param_count);
    else
        param(i)=fixed_params(i);
    end
end

if plot_flag==1
    x_vector=min(x):(max(x)-min(x))/100:max(x);
    figure;
    plot(x,y,'k.',x_vector,f(param(isnan(fixed_params)),x_vector),'r-')
    xlim([min(x) max(x)])
end
end
%% FUNCTION
function cmap = ColorMapDerek(start_col,end_col,invert)
    cmap = [linspace(start_col(1),end_col(1),256)', linspace(start_col(2),end_col(2),256)', linspace(start_col(3),end_col(3),256)'];
    if strcmp(invert,'yes'); cmap = 1-cmap; end
end